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Preface 


Energy  detection  in  the  RF  spectrum  is  the  most  basic  technique  for  signal 
detection.  Typically,  this  requires  establishing  an  energy  detection  threshold  based 
on  a  noise-only  condition  (i.e.,  no  signal  present  in  the  RF  spectrum).  Initially,  the 
area  of  exploration  for  this  report  was  to  examine  the  potential  of  an  automatic 
energy  detection  thresholding  algorithm  based  on  the  RF  measurement.  It  would 
not  require  the  preliminary  RF  spectrum  noise-only  measurements  to  establish  the 
energy  detection  threshold. 

A  series  of  RF  spectrum  measurements  were  collected  at  the  US  Army  Research 
Faboratory  in  2013.  The  data  files  contained  different  segments  of  the  RF  spectrum 
with  various  resolution  bandwidths.  In  addition,  the  data  files  represented  RF 
spectral  conditions  with  and  without  RF  signals.  Initially,  morphological  imaging 
processing  techniques  were  investigated  to  determine  the  potential  for  an  automated 
energy  detection  algorithm. 

The  technique  of  “opening”  was  very  successful  in  determining  an  energy  detection 
threshold  for  the  RF  spectrum  with  signal  and  noise-only  environments.  This  led 
to  the  investigation  of  other  techniques  that  could  be  used  in  automated  energy 
detection  algorithms.  This  exploration  has  resulted  in  a  series  of  techniques  for 
automatic  energy  detection  thresholding  determination.  This  is  the  first  of  5  reports 
that  detail  the  energy  detection  techniques  examined  with  the  recorded  RF  spectrum 
measurements. 14 
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1.  Introduction 


Energy  detection  is  the  simplest  method  of  detecting  a  signal  in  the  frequency 
spectrum.  When  doing  so,  a  comparison  is  made  between  the  frequency  spectral 
component  energy  and  a  detection  threshold  level.  Knowledge  of  the  frequency 
spectrum  is  usually  necessary  to  establish  this  detection  threshold  level.  Various 
methods  can  be  used  to  determine  the  detection  threshold  level,  such  as  establishing 
the  system  noise  statistics  offline  and  setting  the  threshold  for  a  given  probability 
of  detection  versus  probability  of  false  alarm. 

The  goal  of  this  study  was  to  develop  an  algorithm  that  can  analyze  the  frequency 
spectrum  and  establish  a  threshold  value  based  on  the  spectral  components.  The 
automated  processing  employs  techniques  from  image  processing  and  statistical 
analysis.  A  combination  of  techniques  from  these  2  areas  resulted  in  an  algorithm 
for  determining  a  threshold  detection  level.  From  image  processing,  we  applied 
morphological  filtering  to  the  spectral  data  to  establish  a  noise  level.  The  2 
fundamental  operations  in  this  threshold-establishing  algorithm  are  called  “Erode 
and  Dilate”.  This  series  of  operations  is  performed  on  the  spectral  data  until  certain 
statistical  criteria  are  met.  An  offset  was  then  added  to  the  morphological  filtered 
data  to  form  the  threshold  for  energy  detection. 

2.  Data  Collection  and  Statistical  Summary 

The  local  RF  spectrum  was  measured  in  2013  on  the  rooftop  of  building  204  at  the 
US  Army  Research  Laboratory’s  (ARL’s)  Adelphi  location.  An  Agilent  N9342CN 
spectrum  analyzer  and  a  Discone  antenna  were  used  to  collect  RF  spectrum  data. 
This  spectrum  analyzer  was  operated  under  the  control  of  a  Lab  VIEW  software 
program  to  acquire  and  store  data  with  different  resolution  bandwidth  (RBW)  from 
1  kHz  to  1  MHz. 

These  data  files  represent  various  sizes  of  RF  spectral  coverage  from  10  MHz  up 
to  4  GHz.  The  number  of  data  files  for  each  spectral  band  varied.  The  larger  RBW 
files  were  acquired  over  seconds  of  data  acquisition  time  versus  the  small  RBW 
files.  Small  RBW  provides  fine  spectral  resolution,  but  it  impacts  data  acquisition 
time  for  spectral  coverage  and  data  size.  Depending  on  the  RBW  and  spectral 
coverage,  a  data  file  could  require  a  few  hours  of  acquisition. 

Table  1  summarizes  the  data  collection  measurement  files.  The  RF  spectral  bands 
covered  the  spectrum  from  10  MHz  to  4  GHz.  In  general,  various  spectral  bands 
were  measured  with  4  RBW  configuration.  Data  file  size  is  inversely  proportional 
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to  the  RBW.  Data  file  size  is  proportional  to  the  spectral  band  coverage.  The  data 
size  varied  from  approximately  1  KSample  to  4  MSample  data  points  per  file. 


Table  1  Summary  of  the  RF  spectrum  collection 


Spectral  band  coverage 

Number  of  RBW 

measurement 

RBW 

1. 1-1.6  GHz 

4 

1  kHz,  10  kHz,  100  kHz, 

2-3  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

3-4  GHz 

3 

1  MHz 

10  kHz,  100  kHz,  1  MHz 

4-6  GHz 

3 

10  kHz,  100  kHz,  1  MHz 

10  MHz-1  GHz 

4 

1  kHz,  10  kHz,  100  kHz, 

10  MHz-2  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

10  MHz-3  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

10  MHz^l  GHz 

4 

1  MHz 

1  kHz,  10  kHz,  100  kHz, 

100  MHz-1  GHz 

1 

1  MHz 

100  kHz 

3.  Statistical  Processing 

3.1  Statistical  Analysis 

Statistical  analysis  is  the  mathematical  science  dealing  with  the  analysis  or 
interpretation  of  data.  The  data  analyst  uses  a  few  straightforward  statistical 
techniques  as  a  means  of  summarizing  the  collected  data.  These  statistical 
techniques  are  under  the  area  of  descriptive  statistics,  which  is  a  methodology  to 
condense  the  data  in  quantitative  terms. 

In  commercial  prognostics  and  diagnostic  vibrational  monitoring  applications, 
statistical  techniques  that  are  mainly  used  for  alarm  purposes  in  industrial  plants 
are  the  statistical  moments  of  order  2,  3,  and  4.  The  probability  density  function 
(PDF)  of  the  vibrational  time  series  of  a  good  bearing  has  a  Gaussian  distribution 
(also  known  as  a  normal  distribution),  whereas  a  damaged  bearing  results  in  a  non- 
Gaussian  distribution  with  dominant  tails  because  of  a  relative  increase  in  the 
number  of  high  levels  of  acceleration.  These  techniques  can  be  applied  to  the  RF 
spectral  data  with  a  different  interpretation  of  the  results.1 
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3.1.1  Moments 


If  these  moments  are  calculated  about  the  mean,  they  are  called  central  statistical 
moments.  The  first  and  second  moments  are  well  known,  being  the  mean  and  the 
variance,  respectively.  These  are  analogous  to  the  first  and  second  area  moments  of 
inertia  with  the  area  shape  defined  by  the  PDF.  The  third  moment  is  termed 
skewness  and  the  fourth  moment  is  termed  kurtosis.  The  general  equation  for  the 
order  of  moment  is  as  follows: 


N 


Mr 


=  n  2>  - I)P' 

i=  1 


where  p  is  the  order  of  the  moment, 
N  is  the  number  of  data  value, 
i  is  the  index  of  the  data  value,  and 
x  is  the  mean  value  of  the  data  set. 


3.1.2  Mean 

Mean  is  the  most  common  measure  of  a  statistical  distribution.  In  this  case,  mean 
is  the  arithmetic  average  for  a  set  of  measurements. 


i= 1 


3.1.3  Variance 

Variance  is  a  measure  of  the  dispersion  of  a  waveform  about  its  mean — also  called 
the  second  moment  of  the  measurements. 

N 

ff2  =  ~ *)2 ' 
i= 1 


3.1.4  Standard  Deviation 

Standard  deviation  is  a  measure  of  the  variation  of  a  set  of  data  values.  The  standard 
deviation  is  defined  as  the  square  root  of  the  variance  moment. 
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a 


\ 


I 

i= 1 


(Xj  -  X)2 


3.1.5  Kurtosis 

Kurtosis  is  the  fourth  statistical  moment,  normalized  by  the  standard  deviation  to 
the  fourth  power.  It  is  a  measure  of  whether  the  data  are  peaked  or  flat  relative  to  a 
normal  distribution.  The  noise  in  the  RF  spectrum  is  typically  considered  to  have  a 
normal  distribution.  The  normal  distribution  has  a  value  of  3. 

M4 

K  =  —r  . 

a4 

-  *)4 

K~  a* 

N 

k  =  Ysxi  -  ■ 

i= 1 


3.1.6  Maximum 

“Max”  is  the  largest  value  of  a  set  of  numbers. 

y  —  max[x(n)]. 

3.1.7  Minimum 

“Min”  is  the  smallest  value  of  a  set  of  numbers. 

y  =  min[x(n)]. 


3.1.8  Median 

The  statistical  median  is  an  order  statistic  that  gives  the  “middle”  value  of  a  set  of 
samples.  Median  is  the  middle  value  of  a  set  of  data  values  that  divides  the  set  into 
2  groups.  Half  the  groups  exist  below  and  half  exist  above  this  value. 
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3.1.9  Rank  Order  Filter 


Rank  order  filter  is  a  sorting  process  by  which  a  set  of  numbers  is  ordered  from  the 
smallest  to  the  largest  value.  Rank  order  filtering  is  a  nonlinear  filtering  technique 
that  orders  the  contents  of  a  filter  kernel  and  selects  the  sample  indexed  by  rank 
from  the  magnitude  ordered  samples. 

Rank  order  filtering  can  be  summarized  as  follows: 

y  =  Xf=  i  «£*(£)  , 

where  x (q,  l  —  1,  ■  ■  ■ ,  N  is  the  result  of  sorting  the  data  in  ascending  order.  With  this 
definition,  the  max,  min,  and  median  can  be  obtained  from  a  rank  order  filter  as 
follows2: 

Min: 

y  =  £f=  1  UiX(i) 


a,  =  1  for  i  —  1 


0  otherwise 


Max: 


y  =  £f=  1  a£x(0 


cij  =  1  /or  i  =  N 


0  otherwise 


Median: 


y  =  £f=!  a*x(0 


at  —  1  for  i 


N 

~2 


0  otherwise 
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3.1.10  Crest  Factor  (CF) 

Crest  factor  (CF)  is  a  measure  of  a  waveform  showing  the  ratio  of  peak  values  to 
the  effective  value.  In  other  words,  CF  indicates  how  extreme  the  peaks  are  in  a 
waveform. 


Crest  Factor  —  ^  ^peak  _ 

%rms 

Noise  sources  are  characterized  by  their  CF,  which  is  the  peak  to  average  ratio  of 
the  noise.  In  a  technical  bulletin,  XiTRON  reported  CF  values  between  5  and  7  for 
random  noise.3  For  example,  a  5:1  CF  of  the  noise  voltage  is  201og5  =  14  dB.  This 
is  a  measure  of  the  quality  of  the  noise  distributions  and  one  way  to  measure  its 
Gaussian  nature.4  For  the  purpose  of  algorithm  development,  the  CF  equation  was 
modified  as  follows: 


Max 

Crest  Factor  —  - - - . 

Median 


3.1.11  Running  Median 

To  define  the  running  median  filter,  let  {jc  }  be  a  discrete  time  series.  The  running 
median  passes  a  window  over  the  sequence  { jc  }  that  selects,  at  each  instant  m,  an 
odd  number  of  samples  to  comprise  the  observation  vector  x{n).  The  observation 
window  is  typically  symmetric  and  centered  at  n,  resulting  in 

x(n )  =  [x(n  —  Nf),  ■■■  ,x(n),  ,x(n  +  Nf)]7  , 

where  may  range  in  value  over  the  nonnegative  integers,  and  N  —  2Nt  +  1  is 
the  (odd  valued)  window  size.  While  processing  such  noncausal  observation 
vectors  has  traditionally  been  referred  to  as  smoothing,  we  loosen  the  terminology 
somewhat  and  refer  to  the  processing  of  both  causal  and  noncausal  observations  as 
simply  filtering.  The  median  filter  operating  on  the  input  sequence  {x}  produces 
the  output  sequence{y},  where  at  time  index  n, 

y(n)  =  MED  [x(n)] 

=  Median  value  of  [ x(n  —  N- L),  •••  ,x(n),  ,x(n  +  A^)]. 

That  is,  the  samples  in  the  observation  window  are  sorted,  and  the  middle,  or 
median,  value  is  taken  as  the  output.5 


Approved  for  public  release;  distribution  is  unlimited. 

6 


3.2  Statistical  Summary 


There  were  31  different  groupings  of  the  RF  spectrum  data  measurements.  The 
number  of  data  files  under  each  of  the  main  groupings  was  varied.  For  the  purpose 
of  developing  the  algorithm,  only  a  single  data  file  was  selected  from  each  group. 
Each  data  file  was  processed  to  obtain  the  following  characteristics:  RBW,  mean, 
standard  deviation,  median,  max,  min,  kurtosis,  and  CF.  A  summary  of  the  results 
is  in  Table  2. 

The  following  results  were  noted: 

•  The  smaller  the  RBW,  the  lower  the  noise  floor.  The  equation  for  thermal 
noise  power  is  P  —  kTB,  where  B  is  bandwidth.  In  this  case,  the  RBW  of  1 
kHz  has  the  lowest  noise  value.  The  RBW  of  1  MHz  has  the  highest  noise 
value. 

.  Each  10-fold  increase  in  bandwidth  results  in  a  10-dB  increase  in  noise 
power.  This  relationship  is  illustrated  in  this  data  set. 

•  The  calculated  mean  and  median  values  are  very  close  for  a  given  RF 
measurement  configuration. 

•  The  CF  for  noise-only  data  files  was  on  the  order  of  approximately  10  to 
13.  Noise-only  data  files  were  estimated  by  visually  inspecting  the  spectrum 
plot. 

•  The  RF  spectrum  is  not  necessarily  flat  across  the  band. 
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Table  2  Statistical  analysis  of  the  RF  spectrum  measurements 


Filename 

Spectral  Band 

RBW 

Mean 

SD 

Median 

Max 

Min 

Kurtosis  CF 

Ai  r_test_l.  lGHz_1.6GHzb_03_28_14_06_33_22 

1.1-  1.6  GHz 

1000 

-112.694 

3.8444 

-112.3 

-101.2 

-138.4 

3.5571 

11.1 

Ai  r_test_l.  lGHz_l.  6GHza_03_27_14_07_14_05 

1.1-  1.6  GHz 

10000 

-102.12 

3.7226 

-101.8 

-91.69 

-124.5 

3.5134 

10.11 

Ai  r_test_l.  lGHz_l.  6GHzc_03_31_14_06_47_26 

1.1-  1.6  GHz 

100000 

-90.9109 

2.597 

-90.75 

-83.13 

-101.2 

3.0217 

7.62 

Ai  r_test_l.lGHz_l.  6GHzd_04_01_14_06_54_03 

1.1-  1.6  GHz 

1000000 

-80.5634 

3.5741 

-80.745 

-55.36 

-89.8 

14.6861 

25.385 

A  i  r_te  st_2G  H  z_3G  H  za_06_05_14_07_15_58 

2- 3 GHz 

1000 

-111.219 

4.4991 

-111.1 

-77.55 

-139 

6.8918 

33.55 

Ai  r_te  st_2G  H  z_3G  H  z  b_05_29_14_06_28_28 

2- 3  GHz 

10000 

-101.006 

4.1808 

-100.8 

-74.67 

-122 

6.1995 

26.13 

Ai  r_te  st_2G  H  z_3G  H  zc_05_29_14_04_09_27 

2- 3 GHz 

100000 

-89.5666 

3.4154 

-89.68 

-64.27 

-100.3 

13.2282 

25.41 

A  i  r_te  st_2G  H  z_3G  H  zd_05_28_14_00_06_18 

2- 3  GHz 

1000000 

-79.6014 

3.0753 

-79.6 

-58.74 

-88.88 

8.4466 

20.86 

A  i  r_te  st_3G  H  z_4G  H  z  b_06_12_14_06_30_17 

3-4  GHz 

10000 

-100.51 

3.7072 

-100.2 

-89.68 

-123.7 

3.5901 

10.52 

Ai  r_te  st_3G  H  z_4G  H  z  C_06_l  1_14_06_44_40 

3-4  GHz 

100000 

-89.1625 

2.5513 

-89.075 

-80.68 

-102.3 

3.0774 

8.395 

A  i  r_te  st_3G  H  z_4G  H  zd_06_09_14_07_09_5 1 

3-4  GHz 

1000000 

-79.0828 

2.445 

-78.94 

-72.76 

-88.5 

2.9651 

6.18 

A  i  r_te  st_4G  H  z_6G  H  za_04_24_14_07_05_26 

4- 6 GHz 

1000 

-107.559 

4.1793 

-107.2 

-94.31 

-136.3 

3.3696 

12.89 

A  i  r_te  st_4G  H  z_6G  H  z  b_04_28_14_06_23_12 

4- 6  GHz 

10000 

-97.245 

4.0398 

-96.95 

-84.05 

-119.7 

3.3249 

12.9 

A  i  r_te  st_4G  H  z_6G  Hzc_04_29_14_06_48_42 

4- 6  GHz 

100000 

-85.9125 

3.082 

-85.77 

-75.84 

-98.78 

2.9068 

9.93 

A  i  r_te  st_10M  H  z_lG  Hza_04_17_14_07_09_30 

10  MHz  - 1  GHz 

1000 

-110.348 

8.6126 

-112.1 

-28.37 

-140.3 

6.4788 

83.73 

A  i  r_te  st_10M  H  z_lG  Hz  b_04_21_14_07_17_44 

10  MHz  - 1  GHz 

10000 

-100.168 

8.4775 

-102 

-27.95 

-126.7 

6.7602 

74.05 

A  i  r_te  st_10M  H  z_lG  Hzc_04_22_14_06_39_58 

10  MHz  - 1  GHz 

100000 

-87.948 

8.8113 

-90.74 

-27.28 

-103.6 

6.513 

63.46 

Ai  r_te  st_10M  Hz_lG  Hzd_04_23_14_06_36_16 

10  MHz  - 1  GHz 

1000000 

-78.1054 

8.3327 

-80.59 

-27.82 

-89.88 

7.5282 

52.77 

A  i  r_te  st_10M  H  z_2G  H  zd_03_06_14_08_52_04 

10  MHz -2  GHz 

1000 

-111.275 

6.9798 

-112.1 

-27.34 

-138.3 

9.5085 

84.76 

A  i  r_te  st_10M  Hz_2G  H  zC_02_27_14_06_43_55 

10  MHz -2  GHz 

10000 

-100.894 

7.0157 

-101.8 

-26.58 

-125.9 

10.4224 

75.22 

A  i  r_te  st_10M  Hz_2G  H  z  B_02_25_14_06_35_53 

10  MHz -2  GHz 

100000 

-89.242 

6.7487 

-90.59 

-26.09 

-102.4 

12.7961 

64.5 

A  i  r_te  st_10M  H  z_2G  H  z_02_20_14_06_55_36 

10  MHz -2  GHz 

1000000 

-78.8448 

7.3203 

-80.6 

-26.54 

-92.41 

11.3819 

54.06 

Air_test_10MHz_3GHz_03_ll_14_06_32_51 

10  MHz -3  GHz 

1000 

-111.373 

6.1356 

-111.8 

-29.42 

-140.5 

11.0474 

82.38 

Ai  r_te  st_10MHz_3G  Hz  B_03_13_14_06_32_45 

10  MHz -3  GHz 

10000 

-100.787 

6.3591 

-101.4 

-25.36 

-125.9 

11.7254 

76.04 

Air_test_10MHz_3GHzC_03_18_14_06_32_08 

10  MHz -3  GHz 

100000 

-89.1614 

6.0257 

-90.14 

-24.16 

-102.4 

15.3376 

65.98 

Air_test_10MH  z_3G  H  z  d_03_20_14_08_09_07 

10  Mhz  -  3  GHz 

1000000 

-78.6492 

6.7418 

-80.8 

-25.81 

-90.38 

13.0673 

54.99 

Air_test_10MHz_4GHza_04_10_14_07_01_22 

10  MHz -4  GHz 

1000 

-111.158 

5.6813 

-111.4 

-28.22 

-143.1 

10.9573 

83.18 

Air_test_10MHz_4GHzb_04_ll_14_06_10_18 

10  MHz -4  GHz 

10000 

-100.838 

5.5851 

-101.1 

-26.85 

-128 

12.2576 

74.25 

Air_test_10MHz_4GHzc_04_15_14_07_01_09 

10  MHz -4  GHz 

100000 

-89.4431 

5.1076 

-90.03 

-25.05 

-103.5 

18.317 

64.98 

A  i  r_t  e  s  t_10  M  H  z_4G  H  z  d_04_16_14_06_42_27 

10  MHz -4  GHz 

1000000 

-79.0396 

5.5441 

-79.89 

-26.72 

-91.64 

16.9073 

53.17 

A  i  r_t  e  st_100M  Hz_lG  H  z_02_19_14_07_20_46 

100  MHz  - 1  GHz 

100000 

-88.0772 

9.0701 

-90.82 

-23.91 

-102.2 

6.994 

66.91 

4.  Morphological  Image  Processing 

A  description  of  a  technique  for  automatically  estimating  the  noise  floor  spectrum 
was  given  in  a  conference  back  in  1997. 6  This  technique  is  based  on  applying  the 
morphological  binary  image  processing  operators  to  the  RF  spectrum.  The 
technique  works  well  in  both  flat  and  nonflat  noise  floor  spectra.  Morphology 
image  processing  is  a  set  of  nonlinear  operations.  Ready  et  al.  note  that  “humans 
are  good  at  estimating  the  noise  floor  spectrum  by  ‘eyeballing’  a  spectral  plot. 
Intuitively,  we  separate  the  spectral  humps  from  the  noise  floor  spectrum  by 
eliminating  those  parts  of  the  spectrum  shape  that  are  due  to  signals  and  visually 
draw  in  the  noise  floor  spectrum.”6 

Morphology  is  a  broad  set  of  processing  techniques  that  process  images  based  on 
shapes.  Morphological  operations  apply  a  structuring  element  to  an  input  image, 
creating  an  output  image  of  the  same  size.  In  a  morphological  operation,  the  value 
of  each  pixel  in  the  output  image  is  based  on  a  comparison  of  the  corresponding 
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pixel  in  the  input  image  with  its  neighbors.  By  choosing  the  size  and  shape  of  the 
neighborhood,  you  can  construct  a  morphological  operation  that  is  sensitive  to 
specific  shapes  in  the  input  image. 

The  most  basic  morphological  operations  are  dilation  and  erosion.  Dilation  adds 
pixels  to  the  boundaries  of  the  objects  in  an  image,  while  erosion  removes  pixels 
on  object  boundaries.  The  number  of  pixels  added  or  removed  from  the  objects  in 
an  image  depends  on  the  size  and  shape  of  the  structuring  element  used  to  process 
the  image.  In  the  morphological  dilation  and  erosion  operations,  the  state  of  any 
given  pixel  in  the  output  image  is  determined  by  applying  a  rule  to  the 
corresponding  pixel  and  its  neighbors  in  the  input  image.  These  rules  are  known  as 
the  erosion  and  dilation: 

Erosion: 


AQB 


n 

b  E  B 


A 


-ft- 


Dilation: 


A®B 


U  R 
a  E  A  a 


The  opening  of  A  by  B  is  obtained  by  the  dilation  of  A  by  B,  followed  by  erosion  of 
the  resulting  structure  by  B. 

Opening: 


A  o  B  =  (A  ©  B)©B. 

The  application  of  the  opening  technique  on  the  RF  spectrum  is  the  basis  for  noise 
floor  estimation.  The  estimation  of  the  noise  floor  is  used  as  a  reference  level  to  set 
a  threshold  detection  level.7 

In  signal,  statistical,  and  image  processing,  the  minimum  and  maximum  operators 
are  typically  encountered.  Gil  and  Kimmel  note,  “In  mathematical  morphology,  the 
result  of  such  an  operator  is  referred  to  as  the  erosion  (or  dilation)  of  the  signal  with 
a  structuring  element  given  by  a  pulse  of  width  p.”8 

For  the  1-D  case,  this  reduces  to  a  simple  filter  of  just  providing  the  max  or  min 
value  of  a  set  of  values.8 

1-D  max  filter:  Given  a  sequence  x0,  ••• ,  xn_1  and  an  integer  p  >  1,  compute 

V;  =  max  Xj+; 

o<  j<p  J 
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for  =  0,  •••  ,n  —  p  . 


1-D  min  filter:  Given  a  sequence  x0,  ■■■ ,  xn_1  and  an  integer  p  >  1,  compute 

yi  -  min  xi+j 

o  <J<p 

For  =  0 ,n  —  p  . 


5.  Algorithm 

The  algorithm  was  originally  evaluated  and  developed  in  LabVIEW.  This 
development  process  was  then  executed  in  MATLAB.  Appendix  A  is  the  code  used 
to  process  and  generate  the  enclosed  RF  spectrum  signatures  with  the 
corresponding  results  of  the  morphological  processing  and  resultant  threshold 
detection  level.  The  following  is  a  description  of  the  detection  thresholding  level 
generation  based  on  morphological  filter  processing: 

1)  Determine  the  RBW  of  the  spectral  data. 

2)  Determine  some  statistics  on  the  spectral  data  file:  median,  max,  and  CF. 

3)  Perform  a  sliding  running  median  on  data  for  a  window  size  of  5.  The  start 
and  end  of  the  spectral  data  are  somewhat  an  issue.  In  this  case,  the  first  5 
points  used  the  median  value  determined  at  the  first  valid  value.  The  last  5 
points  used  the  median  value  determined  at  the  last  valid  point. 

4)  Perform  an  erosion  operation  on  the  median  filtered  spectral  data  as 
computed  in  step  3.  Starting  from  the  lowest- frequency  component,  create 
the  eroded  spectral  array  with  an  initial  sliding  window  of  size  2.  Basically, 
calculate  the  minimum  value  for  the  given  window,  fill  in  that  position,  and 
shift  the  processing  window  one  frequency  position  to  calculate  the  next 
value.  Continue  on  across  the  entire  median  filtered  spectral  data  array.  For 
the  last  value,  use  the  previous  valid  value  that  is  calculated. 

5)  Perform  a  dilation  operation  on  the  eroded  spectral  array  as  calculated  in 
step  4.  In  this  case,  the  calculations  are  performed  starting  at  the  highest 
frequency  and  continue  to  the  lowest  frequency.  Dilation  is  simply  the 
maximum  value  in  the  windowed  data  set.  The  initial  window  size  is  2.  In 
this  case,  the  maximum  value  is  determined  from  the  sliding  window  data 
and  used  to  create  the  dilated  spectral  array. 

6)  Calculate  statistics  on  the  dilated  spectral  array  from  step  5:  mean,  median, 
max,  and  CF. 
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7)  Repeat  the  morphological  filtering  of  erosion  and  dilation  on  the  spectral 
array  by  increasing  the  window  size  by  1  for  every  pass.  For  example,  on 
the  second  pass,  the  window  size  is  3  for  erosion  and  dilation  operations. 
The  erosion  operation  is  applied  to  the  dilated  spectral  array  from  step  5. 
Repeat  this  process  as  long  as  the  CF  is  greater  than  10  and  the  mean  power 
greater  than  0.01  on  each  pass  of  the  morphological  filtering. 

8)  If  the  CF  is  not  greater  than  10,  then  perform  the  morphological  filtering 
operations  until  the  average  power  is  greater  than  0.01. 

9)  Once  the  convergence  has  been  met,  use  the  resulting  morphological 
filtered  data  array  (i.e.,  the  dilated  spectral  array)  to  form  the  threshold  as 
follows: 

Threshold  — 

Morphological  Filtered  Spectral  Array  +  25  — 

„  _  „  _  ,  (RBW\ 

2.9x10log10(— ) 


6.  Conclusion 


In  general,  the  morphological  filtering  for  establishing  the  thresholding  does  well. 

Appendix  B  displays  the  graphs  of  the  morphological  processing  for  each  of  the 
selected  RF  spectrum  data  files.  The  red  curve  is  the  result  of  the  morphological 
processing  overlaid  on  top  of  the  RF  spectrum  signature. 

Appendix  C  displays  the  results  of  the  threshold  overlaid  on  top  of  the  RF  spectrum 
signature.  The  graphs  are  intended  to  provide  the  reader  with  a  qualitative  sense  of 
the  effectiveness  of  the  automatic  threshold  generation  algorithm.  Actual 
processing  time  and  the  number  of  iterations  through  the  morphological  filtering 
routine  are  highly  dependent  on  the  RBW  and  spectral  coverage  band.  A  small 
RBW  over  a  large  coverage  band  requires  more  time  to  complete  the  process.  The 
threshold  detection  level  obtained  works  well  for  signals  present  or  noise-only  RF 
spectrum  data  files.  The  red  curve  is  the  threshold  when  added  to  the  morphological 
processed  array  as  shown  in  Appendix  B. 

Overall,  a  visual  inspection  of  the  graphs  shows  that  the  algorithm  works  well.  The 
threshold  parameters  can  be  modified  to  obtain  the  desired  level  of  detection  versus 
the  false  alarm.  The  algorithm  worked  well  for  the  RF  spectrum  where  there  are 
signals  present  as  well  as  for  the  RF  spectrum  that  was  judged  to  be  noise  only. 
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Appendix  A.  MATLAB  Code 
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function  MorphologicalFilterR0F9report  () 


%Make  selection  for  data  file  to  process 
%SELECT  =  1  to  31 

SELECT  =  15; 


dir{l}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHza\ ' ; 

filename{l}  =  ' Air_test_l . lGHz_l . 6GHza_03_27_14_07_14_05 ' ; 

dir{2}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHzb\ '  ; 

filename{2}  =  ' Air_test_l . lGHz_l . 6GHzb_03_28_14_06_33_22 ' ; 

dir{3}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHzc\ ' ; 

filename{3}  =  'Air_test_l . lGHz_l . 6GHzc_03_31_14_0 6_47_2 6 ' ; 

dir{4}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_l . lGHz_l . 6GHzd\ '  ; 

filename{4}  =  'Air_test_l . lGHz_l . 6GHzd_04_01_14_06_54_03 ' ; 

dir{5}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHza\ ' ; 

filename{5}  =  'Air_test_2GHz_3GHza_06_05_14_07_15_58 ' ; 

dir{6}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHzb\ '  ; 

filename{6}  =  ' Air_test_2GHz_3GHzb_0 5_2 9_1 4_0 6_2 8_2 8 ' ; 

dir{7}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHzc\ ' ; 

filename{7}  =  ' Air_test_2GHz_3GHzc_05_29_14_04_09_27 ' ; 

dir{8}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_2GHz_3GHzd\ '  ; 

filename{8}  =  ' Air_test_2GHz_3GHzd_05_28_14_00_06_18 ' ; 

dir{9}  =  'K:\CognitiveRadar\spectrum  monitoring\data\building204- 
4c085\Air_test_3GHz_4GHzb\ ' ; 

filename{9}  =  'Air_test_3GHz_4GHzb_06_12_14_06_30_17 ' ; 
dir{10}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_3GHz_4GHzC\ ' ; 
filename{10}  =  ' Air_test_3GHz_4GHzC_06_ll_14_06_44_40 ' ; 

dir{ll}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_3GHz_4GHzd\ ' ; 
filename{ll}  =  ' Air_test_3GHz_4GHzd_06_09_14_07_09_51 ' ; 
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dir{12}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHza\ ' ; 
filename{12}  =  ' Air_test_4GHz_6GHza_04_24_14_07_05_26 ' ; 

dir{13}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHzb\ ' ; 
filename{13}  =  ' Air_test_4GHz_6GHzb_04_28_14_06_23_12 '  ; 

dir{14}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHzc\ ' ; 
filename{14}  =  ' Air_test_4GHz_6GHzc_04_29_14_06_48_42 ' ; 

dir{15}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c085\Air_test_10MHz_lGHza\ ' 
filename{15}  =  ' Air_test_10MHz_lGHza_04_17_14_07_09_30 ' ; 

dir{16}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c0 8 5\Air_test_10MHz_lGHzb\ ' 
filename{16}  =  ' Air_test_10MHz_lGHzb_04_21_14_07_17_44 ' ; 

dir{17}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_lGHzc\ ' 
filename{17}  =  'Air_test_10MHz_lGHzc_04_22_14_06_39_58 ' ; 

dir{18}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_lGHzd\ ' 
filename{18}  =  ' Air_test_10MHz_lGHzd_04_23_14_06_36_16 ' ; 

dir{19}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHz\ ' ; 
filename{19}  =  'Air_test_10MHz_2GHz_02_20_14_06_55_36 '  ; 

dir{20}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHzB\ ' 
filename{20}  =  ' Air_test_10MHz_2GHzB_02_25_14_06_35_53 ' ; 

dir{21}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c0 8 5\Air_test_10MHz_2GHzC\ ' 
filename{21}  =  ' Air_test_10MHz_2GHzC_02_27_14_06_43_55 ' ; 

dir{22}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c085\Air_test_10MHz_2GHzd\ ' 
filename{22}  =  ' Air_test_10MHz_2GHzd_03_06_14_08_52_04 ' ; 

dir{23}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHz\ ' ; 
filename{23}  =  ' Air_test_10MHz_3GHz_03_ll_14_06_32_51 ' ; 

dir{24}  =  'K:\CognitiveRadar\spectrum 

monitor ing\data\building204-4c0 8 5\Air_test_10MHz_3GHzB\ ' 
filename{24}  =  'Air_test_10MHz_3GHzB_03_13_14_06_32_45 ' ; 

dir{25}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHzC\ ' 
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filename { 25 } 


=  ' Air_test_10MHz_3GHzC_03_18_14_06_32_08 ' ; 
dir{26}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHzd\ ' ; 
filename{26}  =  ' Air_test_10MHz_3GHzd_03_20_14_08_09_07 ' ; 

dir{27}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHza\ ' ; 
filename{27}  =  ' Air_test_10MHz_4GHza_04_10_14_07_01_22 ' ; 

dir{28}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzb\ ' ; 
filename{28}  =  ' Air_test_10MHz_4GHzb_04_ll_14_06_10_18 ' ; 

dir{29}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzc\ ' ; 
filename{29}  =  ' Air_test_10MHz_4GHzc_04_15_14_07_01_09 ' ; 

dir{30}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzd\ ' ; 
filename{30}  =  ' Air_test_10MHz_4GHzd_04_16_14_06_42_27 ' ; 

dir{31}  =  'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_100MHz_lGHz\ ' ; 
filename{31}  =  ' Air_test_100MHz_lGHz_02_19_14_07_20_46 ' ; 

str_meta  = 

sprint f ( ' %s%s . mspecrawdata ' , dir { SELECT } , filename { SELECT } ) 
str_data  = 

sprint f ( ' %s%s . specrawdata ' , dir { SELECT } , filename { SELECT } ) ; 
%  Get  metadata  on  selected  data  file 


fid__meta  =  f open  ( str_meta)  ; 

META  =  textscan (fid_meta, '%s'); 

ave  =  META { 1 } { 1 }  ; 

ref  =  META{ 1 } {2 }  ; 

attn  =  META{ 1 } {3}; 

rbw  =  META{ 1 } {4}; 

%  Read  in  data  file 

fid_data  =  f open ( str_data) ; 

g=0; 

f=0; 

a=0 ; 

while (g==0 ) 

ft=fgetl (fid_data)  ; 
f = [ f , str2num ( ft )  ]  ; 
at=fgetl (fid_data)  ; 
a= [ a, str2num (at ) ] ; 
g=feof (fid_data)  ; 

end 
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f=f (2 : end) ; 
a=a (2 : end) ; 

%  Determine  resolution  bandwidth 

Bandwidth  =  f(3)  -  f(2); 

%  Plot  data  file 

figure ( 1 ) ; 
plot (f /le6, a) ; 
axis  tight; 
grid; 

xlabel ( 'Frequency  (MHz)  ' )  ; 

ylabel ( 'Power  (dBm)  ' )  ; 

title  (strcat (rbw,  ' , ave,  '  —  ' , ref ,  '  - 

' , attn) , ' Interpreter ' , ’ none ' ) ; 


%  Determine  some  statistics  values  for  data 
M  =  mean (a) ; 

Med  =  median (a); 

S  =  std (a) ; 

Max  =  max (a) ; 

Min  =  min  (a) ; 

Kurt  =  kurtosis  (a) ; 

Range  =  abs (Min  -  Max) ; 

Number  =  floor (Range ) ; 

Bins  =  Number  *  2; 

disp ( [M  S  Med  Max  Min  Kurt] )  ; 


%  Initialize  parameters 

Inputlndex  =  length (a) ; 
Erosionlndex  =  1; 
InputArray  =  a; 
ErosionArray  =  a; 
DilationArray  =  a; 
Medsize  =  5; 

Medcount  =  fix (Medsize) ; 
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Medstart  =  Medcount  +  1; 

%  Process  data  file  with  a  running  median 

for  m  =  Medstart : Inputlndex-Medcount 

InputArray (m)  =  median (a (m-Medcount : m+Medcount ))  ; 

end 

for  m  =  1: Medcount 

InputArray (m)  =  InputArray (Medstart ) ; 

end 

for  m  =  Inputlndex-Medcount : Inputlndex 

InputArray (m)  =  InputArray ( Inputlndex-Medcount ) ; 

end 

%  Plot  processed  median  of  data  file 

figure (2 ) ; 

plot (InputArray) ; 

PreviousPower  =  sum (a) /Inputlndex; 
display  (PreviousPower) ; 

%  Initialize  parameters 

j=  i; 

DiffPower  =  1; 

DilationMedian  =  Med; 

DilationMax  =  Max; 

DilationCF  =  DilationMax-DilationMedian; 


if  DilationCF  >  10  &&  j  ==  1 
display ( ' CF  &  Power'); 

%  Perform  Erosion  operation  on  data  file 

while  DilationCF  >  10 

for  k  =  1 : Inputlndex- j 

ErosionArray (k)  =  min ( InputArray (k : k+j )) ; 

end 


for  k  =  Inputlndex- j : Inputlndex 
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ErosionArray (k)  =  InputArray ( Input Index- j ) ; 


end 

%  Plot  results  of  Erosion  operation 

figure  (5) ; 

plot (ErosionArray) ; 

%  Perform  Dilation  operation  on  data  file 
DilationArray  =  ErosionArray; 
for  k  =  Inputlndex : j+1 

DilationArray (k)  =  max (ErosionArray (k-j : k) ) ; 

end 


for  k  =  j : 1 

DilationArray (k)  =  ErosionArray ( j+1 ) ; 

end 

%  Plot  results  of  Dilation  operation 

InputArray  =  DilationArray; 
figure  (6) ; 

plot (DilationArray) ; 

DilationPower  =  sum (DilationArray) /Inputlndex; 
DilationMedian  =  median (DilationArray) ; 

DilationMax  =  max (DilationArray)  ; 

DiffPower  =  abs  (PreviousPower  -  DilationPower) ; 
DilationCF  =  DilationMax-DilationMedian; 

display  ([j  DilationPower  DilationMedian  DilationMax 
DilationCF  DiffPower] ) ; 

PreviousPower  =  DilationPower; 

j=  j+i; 

end 

%  Check  performance  of  the  Erosion  &  Dilation  to  see  if 
operations 

%  resulted  have  converged  to  accept  criteria  (If  power 
difference  in 
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%  the  RF  spectrum  is  less  than  0.1  between  iterations  previous 
and 

%  current  processing  of  RF  spectrum 
while  DiffPower  >  0.1 


display ( ' Power ' ) ; 


for  k  =  1 : Inputlndex- j 

ErosionArray (k)  =  min ( InputArray (k : k+ j ) )  ; 


end 

for  k  =  Inputlndex- j : Inputlndex 

ErosionArray (k)  =  min ( InputArray (k : Inputlndex) )  ; 

end 


figure ( 5 ) ; 

plot (ErosionArray) ; 


DilationArray  =  ErosionArray; 


for  k  =  Inputlndex : j+1 

DilationArray (k)  =  max (ErosionArray (k-j : k) )  ; 


end 

for  k  =  j : 1 

DilationArray (k)  =  max (ErosionArray (k : 1 )) ; 

end 

InputArray  =  DilationArray; 
figure (6) ; 

plot  (DilationArray) ; 

DilationPower  =  sum (DilationArray) /Inputlndex; 
DilationMedian  =  median (DilationArray) ; 

DilationMax  =  max (DilationArray) ; 

DiffPower  =  abs  (PreviousPower  -  DilationPower) ; 
DilationCF  =  DilationMax-DilationMedian; 
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display  ([j  DilationPower  DilationMedian  DilationMax 
DilationCF  DiffPower] ) ; 

PreviousPower  =  DilationPower; 

j=  j+i; 

end 

else 


while  DiffPower  >  0.01 
display ( ' Power ' ) ; 


for  k  =  1 :  Inputlndex- j 

ErosionArray (k)  =  min ( InputArray (k : k+ j ) ) ; 


end 

for  k  =  Inputlndex- j : Inputlndex 

ErosionArray (k)  =  min ( InputArray (k : Inputlndex) ) 

end 


figure  (5) ; 

plot (ErosionArray) ; 


DilationArray  =  ErosionArray; 


for  k  =  Inputlndex : j+1 

DilationArray (k)  =  max (ErosionArray (k-j : k) ) ; 


end 

for  k  =  j : 1 

DilationArray (k)  =  max (ErosionArray (k : 1 )) ; 

end 

InputArray  =  DilationArray; 
figure (6) ; 

plot  (DilationArray) ; 

DilationPower  =  sum (DilationArray) /Inputlndex; 
DilationMedian  =  median (DilationArray) ; 

DilationMax  =  max (DilationArray); 
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DiffPower  =  abs  (PreviousPower  -  DilationPower) ; 
DilationCF  =  DilationMax-DilationMedian; 

display  ( [ j  DilationPower  DilationMedian  DilationMax 
DilationCF  DiffPower] ) ; 

PreviousPower  =  DilationPower; 

j=  j+i; 

end 


end 

%  Plot  original  RF  spectrum  and  threshold 

figure  ( 3 ) ; 

plot (f /le6, a) ; 
axis  tight; 

grid; 

xlabel ( 'Frequency  (MHz)  ' )  ; 
ylabel ( 'Power  (dBm) ' ) ; 

title  (strcat (filename (SELECT) ,  '  -  -  ' , rbw,  ' 

Hz'),'  Interpreter ' ,  ' none ' ) ; 
hold  on 

%  Determine  threshold  array 

Threshold  =  DilationArray+25- (2 . 9*logl0 (Bandwidth)); 


plot  (f/le6, Threshold,  :r'); 


hold  off; 

%  Plot  original  RF  spectrum  and  final  results  of  Erosion  & 
Dilation 
%  processing 

figure ( 10 ) ; 

plot  (f /le6, a) ; 
axis  tight; 

grid; 

xlabel ( 'Frequency  (MHz)  ' )  ; 
ylabel ( 'Power  (dBm) ' ) ; 
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rbw 


title  (strcat (filename (SELECT) ,  ' 
Hz'), 'Interpreter ' , ' none ' ) ; 
hold  on 


plot (f/le6, DilationArray , ' r ' ) ; 


hold  off; 

%  Save  results  to  data  file 

lab  =  num2str (SELECT)  ; 
labl  =  strcat (' R0F2-Figure  '  , lab,  '  . 
saveas (gcf ,  labl); 

end 
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jpg' ) ; 


Intentionally  left  blank. 
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Appendix  B.  Graphs  of  Morphological  Processed  RF  Spectrum 

Files 
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Appendix  C.  Graphs  of  RF  Spectrum  Files  Calculated  Detection 

Threshold 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


ARL 

US  Army  Research  Laboratory 

CF 

crest  factor 

GHz 

gigahertz 

kHz 

kilohertz 

MHz 

megahertz 

PDF 

probability  density  function 

RBW 

resolution  bandwidth 

RF 

radio  frequency 
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